##################################################################
# File:      replication_code_main_analysis.R
# Purporse:  Replication of the main text figures and tables for 
#                 Journal of Conflict Resolution
# Author:    Nadiya Kostyuk and Erik Gartzke
# Last Edit: 01/05/2023
##################################################################

rm(list=ls())

## Install & load packages (all at once)
list.of.packages <- c('knitr', 'countrycode', 'dplyr','foreign','glmmTMB','lme4',
                      'stargazer', 'doBy', 'ggplot2', 'sjPlot')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}
lapply(list.of.packages, require, character.only = TRUE)

# setting directory: 
setwd("")


# loading data: 
load('main_data.Rdata')

# 1. covariates, no interactions
mod1 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA)
                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                        +CINC_A_sc+INV_CAP_DISTANCE_sc
                        +NUCLEAR
                        +(1+TYPE2||DYAD),
                        data=data,
                        family=binomial,
                        control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod1)


# 2. adding interactions with covariates (i.e., internet users per capita for the attacker and the target)
mod2 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc)
                          +CINC_A_sc+INV_CAP_DISTANCE_sc
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod2)


# 3. adding interactions with covariates (i.e., internet users per capita for the attacker and the target, cinc score for the attacker)
mod3 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc)+INV_CAP_DISTANCE_sc
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod3)


# 4. adding interactions with covariates (i.e., internet users per capita for the attacker and the target, cinc score for the attacker, distance between capitals)
mod4 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc+INV_CAP_DISTANCE_sc)
                          +NUCLEAR
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod4)


# 5. adding interactions with all covariates 
mod5 <- lme4::glmer(OUTCOME~TYPE*(CYBER_LAG_NA+CONVENTIONAL_LAG_NA
                                        +LOG_INT_USERS_A_sc+LOG_INT_USERS_B_sc
                                        +CINC_A_sc+INV_CAP_DISTANCE_sc
                                        +NUCLEAR)
                          +(1+TYPE2||DYAD),
                          data=data,
                          family=binomial,
                          control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e5))
)
summary(mod5)



##############################
# Estimating Contrasts
##############################
# Creating Table 1

###################
# mod1
    
L2 <- matrix(0, nrow=11, ncol = length(fixef(mod1)))
L2
colnames(L2) = names(fixef(mod1))
rownames(L2) = c('P_cyber_base', 'P_military_base',
                 'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                 'OR_military_mil(lag)', 'OR_cyber_mil(lag)',
                 'OR_conflict_intus(a)','OR_conflict_intus(b)',
                 'OR_conflict_cinc(a)', 'OR_conflict_distance',
                 'OR_conflict_nuclear(b)') # OR (odds ratio)
L2[1,c('(Intercept)', 'TYPE')] <- c(1,1)
L2[2,c('(Intercept)', 'TYPE')] <- c(1,-1)
L2[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1)
L2[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
L2[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
L2[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
L2[7,c('LOG_INT_USERS_A_sc')] <- 1
L2[8,c('LOG_INT_USERS_B_sc')] <- 1
L2[9,c('CINC_A_sc')] <- 1
L2[10,c('INV_CAP_DISTANCE_sc')] <- 1
L2[11,c('NUCLEAR')] <- 1
L2
res1 <- doBy::esticon(mod1, L=L2)
res1 

res1 = res1 %>%
  dplyr::mutate(starz=case_when(
    p.value>=.1~'',
    p.value>=0.05&p.value<0.1~'^',
    p.value>=0.01&p.value<0.05~'*',
    p.value>=0.001&p.value<0.01~'**',
    p.value<0.001~'***'
  ))
res1

m=1.96
res1_1 <- with(res1[1:2,],
               sprintf('%4.2f%s (%4.2f, %4.2f)',
                       plogis(estimate),
                       starz, 
                       plogis(estimate-m*std.error),
                       plogis(estimate+m*std.error)
               ) 
           
)
res1_1
names(res1_1) = rownames(L2[1:2,])
res1_1
res1_2 <- with(res1[3:11,],
               sprintf('%4.2f%s (%4.2f, %4.2f)',
                       exp(estimate),
                       starz,
                       exp(estimate-m*std.error),
                       exp(estimate+m*std.error)
               ) 
               
)
res1_2
names(res1_2) = rownames(L2[3:11,])
res1_2
res1_2 = t(res1_2)
res1_2 = t(res1_2)
res1_2


aic = round(AIC(mod1), 1)
nobz = nobs(mod1)
res1_2 = rbind(res1_2, aic, nobz)
res1_2


########################
# mod2

L <- matrix(0, nrow=15, ncol = length(summary(mod2)$coefficients[,1]))
L
colnames(L) = names(summary(mod2)$coefficients[,1])
rownames(L) = c('P_cyber_base', 'P_military_base',
                'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                'OR_military_mil(lag)', 'OR_cyber_mil(lag)',
                'OR_conflict_intus(a)','OR_conflict_intus(b)',
                'OR_conflict_cinc(a)', 'OR_conflict_distance',
                'OR_conflict_nuclear(b)',
                'OR_cyber_intus(a)','OR_military_intus(a)',
                'OR_cyber_intus(b)','OR_military_intus(b)') 
L[1,c('(Intercept)', 'TYPE')] <- c(1,1) 
L[2,c('(Intercept)', 'TYPE')] <- c(1,-1) 
L[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1) 
L[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
L[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
L[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
L[7,c('LOG_INT_USERS_A_sc')] <- 1
L[8,c('LOG_INT_USERS_B_sc')] <- 1
L[9,c('CINC_A_sc')] <- 1
L[10,c('INV_CAP_DISTANCE_sc')] <- 1 
L[11,c('NUCLEAR')] <- 1
L[12,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,1) 
L[13,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,-1)
L[14,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,1) 
L[15,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,-1)
res <- doBy::esticon(mod2, L=L)
res

res = res %>%
  dplyr::mutate(starz=case_when(
    p.value>=.1~'',
    p.value>=0.05&p.value<0.1~'^',
    p.value>=0.01&p.value<0.05~'*',
    p.value>=0.001&p.value<0.01~'**',
    p.value<0.001~'***'
  ))
res

m=1.96
res2_1 <- with(res[1:2,],
               sprintf('%4.2f (%4.2f, %4.2f)',
                       plogis(estimate),
                       plogis(estimate-m*std.error),
                       plogis(estimate+m*std.error)
               ) 
               
) 
res2_1
names(res2_1) = rownames(L[1:2,])
res2_1

res2_2 <- with(res[3:15,],
               sprintf('%4.2f%s(%4.2f;%4.2f)',
                       exp(estimate),
                       starz,
                       exp(estimate-m*std.error),
                       exp(estimate+m*std.error)
               ) 
               
) 
res2_2
names(res2_2) = rownames(L[3:15,])
res2_2

aic = round(AIC(mod2), 1)
nobz = nobs(mod2)
res2_2 = c(res2_2, aic = aic, nobz = nobz)



####################################
# mod3

L <- matrix(0, nrow=17, ncol = length(summary(mod3)$coefficients[,1]))
L
colnames(L) = names(summary(mod3)$coefficients[,1])
rownames(L) = c('P_cyber_base', 'P_military_base',
                'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                'OR_military_mil(lag)', 'OR_cyber_mil(lag)',
                'OR_conflict_intus(a)','OR_conflict_intus(b)',
                'OR_conflict_cinc(a)', 'OR_conflict_distance',
                'OR_conflict_nuclear(b)',
                'OR_cyber_intus(a)','OR_military_intus(a)',
                'OR_cyber_intus(b)','OR_military_intus(b)',
                'OR_cyber_cinc(a)', 'OR_military_cinc(a)') 
L[1,c('(Intercept)', 'TYPE')] <- c(1,1) 
L[2,c('(Intercept)', 'TYPE')] <- c(1,-1) 
L[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1) 
L[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
L[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
L[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
L[7,c('LOG_INT_USERS_A_sc')] <- 1
L[8,c('LOG_INT_USERS_B_sc')] <- 1
L[9,c('CINC_A_sc')] <- 1
L[10,c('INV_CAP_DISTANCE_sc')] <- 1 
L[11,c('NUCLEAR')] <- 1
L[12,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,1) 
L[13,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,-1)
L[14,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,1) 
L[15,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,-1)
L[16,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,1) 
L[17,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,-1)
L
res <- doBy::esticon(mod3, L=L)
res

res = res %>%
  dplyr::mutate(starz=case_when(
    p.value>=.1~'',
    p.value>=0.05&p.value<0.1~'^',
    p.value>=0.01&p.value<0.05~'*',
    p.value>=0.001&p.value<0.01~'**',
    p.value<0.001~'***'
  ))
res

m=1.96
res3_1 <- with(res[1:2,],
               sprintf('%4.2f (%4.2f, %4.2f)',
                       plogis(estimate),
                       plogis(estimate-m*std.error),
                       plogis(estimate+m*std.error)
               ) # replace prob with 2 decimal places
               
) 
res3_1
names(res3_1) = rownames(L[1:2,])
res3_1

res3_2 <- with(res[3:17,],
               sprintf('%4.2f%s(%4.2f;%4.2f)',
                       exp(estimate),
                       starz,
                       exp(estimate-m*std.error),
                       exp(estimate+m*std.error)
               ) 
               
) 
res3_2
names(res3_2) = rownames(L[3:17,])
res3_2


aic = round(AIC(mod3), 1)
nobz = nobs(mod3)
res3_2 = c(res3_2, aic = aic, nobz = nobz)
      
      


####################################
# mod4

L <- matrix(0, nrow=19, ncol = length(summary(mod4)$coefficients[,1]))
L
colnames(L) = names(summary(mod4)$coefficients[,1])
rownames(L) = c('P_cyber_base', 'P_military_base',
                'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                'OR_military_mil(lag)', 'OR_cyber_mil(lag)',
                'OR_conflict_intus(a)','OR_conflict_intus(b)',
                'OR_conflict_cinc(a)', 'OR_conflict_distance',
                'OR_conflict_nuclear(b)',
                'OR_cyber_intus(a)','OR_military_intus(a)',
                'OR_cyber_intus(b)','OR_military_intus(b)',
                'OR_cyber_cinc(a)', 'OR_military_cinc(a)', 
                'OR_cyber_distance', 'OR_military_distance'
) 
L[1,c('(Intercept)', 'TYPE')] <- c(1,1) 
L[2,c('(Intercept)', 'TYPE')] <- c(1,-1) 
L[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1) 
L[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
L[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
L[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
L[7,c('LOG_INT_USERS_A_sc')] <- 1
L[8,c('LOG_INT_USERS_B_sc')] <- 1
L[9,c('CINC_A_sc')] <- 1
L[10,c('INV_CAP_DISTANCE_sc')] <- 1
L[11,c('NUCLEAR')] <- 1
L[12,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,1) 
L[13,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,-1)
L[14,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,1) 
L[15,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,-1)
L[16,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,1) 
L[17,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,-1)
L[18,c('INV_CAP_DISTANCE_sc', 'TYPE:INV_CAP_DISTANCE_sc')] <- c(1,1) 
L[19,c('INV_CAP_DISTANCE_sc', 'TYPE:INV_CAP_DISTANCE_sc')] <- c(1,-1)
L
res <- doBy::esticon(mod4, L=L)
res

res = res %>%
  dplyr::mutate(starz=case_when(
    p.value>=.1~'',
    p.value>=0.05&p.value<0.1~'^',
    p.value>=0.01&p.value<0.05~'*',
    p.value>=0.001&p.value<0.01~'**',
    p.value<0.001~'***'
  ))
res

m=1.96
res4_1 <- with(res[1:2,],
               sprintf('%4.2f (%4.2f, %4.2f)',
                       plogis(estimate),
                       plogis(estimate-m*std.error),
                       plogis(estimate+m*std.error)
               ) 
               
) 
res4_1
names(res4_1) = rownames(L[1:2,])
res4_1

res4_2 <- with(res[3:19,],
               sprintf('%4.2f%s(%4.2f;%4.2f)',
                       exp(estimate),
                       starz,
                       exp(estimate-m*std.error),
                       exp(estimate+m*std.error)
               ) 
               
) 
res4_2
names(res4_2) = rownames(L[3:19,])
res4_2


# adding AIC:
aic = round(AIC(mod4), 1)
nobz = nobs(mod_list[[14]])
res4_2 = rbind(res3_2, aic, nobz)
    
 
####################################
# mod5

L <- matrix(0, nrow=21, ncol = length(summary(mod5)$coefficients[,1]))
L
colnames(L) = names(summary(mod5)$coefficients[,1])
rownames(L) = c('P_cyber_base', 'P_military_base',
                'OR_cyber_cyb(lag)', 'OR_military_cyb(lag)',
                'OR_military_mil(lag)', 'OR_cyber_mil(lag)',
                'OR_conflict_intus(a)','OR_conflict_intus(b)',
                'OR_conflict_cinc(a)', 'OR_conflict_distance',
                'OR_conflict_nuclear(b)',
                'OR_cyber_intus(a)','OR_military_intus(a)',
                'OR_cyber_intus(b)','OR_military_intus(b)',
                'OR_cyber_cinc(a)', 'OR_military_cinc(a)', 
                'OR_cyber_distance', 'OR_military_distance',
                'OR_cyber_nuclear(b)', 'OR_military_nuclear(b)'
) 
L[1,c('(Intercept)', 'TYPE')] <- c(1,1) 
L[2,c('(Intercept)', 'TYPE')] <- c(1,-1) 
L[3,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,1) 
L[4,c('CYBER_LAG_NA', 'TYPE:CYBER_LAG_NA')] <- c(1,-1)
L[5,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,-1)
L[6,c('CONVENTIONAL_LAG_NA', 'TYPE:CONVENTIONAL_LAG_NA')] <- c(1,1)
L[7,c('LOG_INT_USERS_A_sc')] <- 1
L[8,c('LOG_INT_USERS_B_sc')] <- 1
L[9,c('CINC_A_sc')] <- 1
L[10,c('INV_CAP_DISTANCE_sc')] <- 1
L[11,c('NUCLEAR')] <- 1
L[12,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,1) 
L[13,c('LOG_INT_USERS_A_sc', 'TYPE:LOG_INT_USERS_A_sc')] <- c(1,-1)
L[14,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,1) 
L[15,c('LOG_INT_USERS_B_sc', 'TYPE:LOG_INT_USERS_B_sc')] <- c(1,-1)
L[16,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,1) 
L[17,c('CINC_A_sc', 'TYPE:CINC_A_sc')] <- c(1,-1)
L[18,c('INV_CAP_DISTANCE_sc', 'TYPE:INV_CAP_DISTANCE_sc')] <- c(1,1) 
L[19,c('INV_CAP_DISTANCE_sc', 'TYPE:INV_CAP_DISTANCE_sc')] <- c(1,-1)
L[20,c('NUCLEAR', 'TYPE:NUCLEAR')] <- c(1,1) 
L[21,c('NUCLEAR', 'TYPE:NUCLEAR')] <- c(1,-1)
L
res <- doBy::esticon(mod5, L=L)
res

res = res %>%
  dplyr::mutate(starz=case_when(
    p.value>=.1~'',
    p.value>=0.05&p.value<0.1~'^',
    p.value>=0.01&p.value<0.05~'*',
    p.value>=0.001&p.value<0.01~'**',
    p.value<0.001~'***'
  ))
res

m=1.96
res5_1 <- with(res[1:2,],
               sprintf('%4.2f (%4.2f, %4.2f)',
                       plogis(estimate),
                       plogis(estimate-m*std.error),
                       plogis(estimate+m*std.error)
               ) 
               
) 
res5_1
names(res5_1) = rownames(L[1:2,])
res5_1

res5_2 <- with(res[3:21,],
               sprintf('%4.2f%s(%4.2f;%4.2f)',
                       exp(estimate),
                       starz,
                       exp(estimate-m*std.error),
                       exp(estimate+m*std.error)
               ) 
               
) 
res5_2
names(res3_2) = rownames(L[3:21,])
res5_2


# adding AIC:
aic = round(AIC(mod5), 1)
nobz = nobs(mod5)
res5_2 = rbind(res5_2, aic, nobz)



#######################################
# combining all models into one table:

d1 = tibble(effect = rownames(res1_2), model_1 = res1_2[,1])
d2 = tibble(effect = rownames(res2_2), model_2 = res2_2[,1])
d3 = tibble(effect = rownames(res3_2), model_3 = res3_2[,1])
d4 = tibble(effect = rownames(res4_2), model_4 = res4_2[,1])
d5 = tibble(effect = rownames(res5_2), model_5 = res5_2[,1])


table = d1 %>%
  full_join(d2, by = 'effect') %>%
  full_join(d3, by = 'effect') %>%
  full_join(d4, by = 'effect') %>%
  full_join(d5, by = 'effect')

writeLines(capture.output(xtable(table)), paste0('main_all.tex'))


#################################################################
# Note: to create Table 2, please use the same approach as above 
# specifically (mod 2 and mod 3)
# but load the results for the espionage and disrupt-degrad models
##################################################################


############################################
# Interaction Plots
# Figure 1
############################################

# Note: using the main model here (mod3)

################
# loading data:
load('main_data.RData')

set_theme(
  base = theme_classic(), axis.linecolor = "grey50",
  axis.textcolor = "black") 
plot_model(mod2, type = "pred", terms = c("LOG_INT_USERS_A_sc[all]", "TYPE"),
           legend.title = '')+
  ggtitle("") +
  scale_color_grey(labels = c("Military", "Cyber"))+
  scale_fill_grey()+
  xlab("\n Attacker's Internet Users per Capita (log, sc)") + 
  ylab("Predicted probability of engaging in conflict \n")

plot_model(mod2, type = "pred", terms = c("LOG_INT_USERS_B_sc[all]", "TYPE"),
           #colors='gs', 
           legend.title = '')+
  ggtitle("") +
  scale_color_grey(labels = c("Military", "Cyber"))+
  scale_fill_grey()+
  xlab("\n Target's Internet Users per Capita (log, sc)") + 
  ylab("Predicted probability of engaging in conflict \n")












